#pragma once

#include <string>
#include <vector>
#include <iostream>
#include <algorithm>

#include "../math/vec3.h"
#include "../math/io.h"

namespace roy {

    namespace mesh {

        template <typename T>
            class tmesh
        {
        public:
            tmesh(void)
            {
                this->Init();
            }

            tmesh(const tmesh& tm)
            {
                this->_v_2_e          = tm._v_2_e;
                this->_v_2_f          = tm._v_2_f;
                this->_e_2_v          = tm._e_2_v;
                this->_e_2_f          = tm._e_2_f;
                this->_f_2_v          = tm._f_2_v;
                this->_f_2_e          = tm._f_2_e;

                this->_v_cnt          = tm._v_cnt;
                this->_e_cnt          = tm._e_cnt;
                this->_f_cnt          = tm._f_cnt;

                this->_v_pos          = tm._v_pos;
                this->_v_n            = tm._v_n;
                this->_e_len          = tm._e_len;
                this->_e_len_average  = tm._e_len_average;
                this->_f_n            = tm._f_n;
                this->_f_center       = tm._f_center;
                this->_f_size         = tm._f_size;

                this->_center                   = tm._center;
                this->_bounding_box[0]          = tm._bounding_box[0];
                this->_bounding_box[1]          = tm._bounding_box[1];
                this->_bounding_box_diagonal    = tm._bounding_box_diagonal;
                this->_bounding_sphere_radius   = tm._bounding_sphere_radius;
            }

            virtual ~tmesh(void)
            {
            }

            tmesh& operator = (const tmesh& tm)
            {
                this->_v_2_e          = tm._v_2_e;
                this->_v_2_f          = tm._v_2_f;
                this->_e_2_v          = tm._e_2_v;
                this->_e_2_f          = tm._e_2_f;
                this->_f_2_v          = tm._f_2_v;
                this->_f_2_e          = tm._f_2_e;

                this->_v_cnt          = tm._v_cnt;
                this->_e_cnt          = tm._e_cnt;
                this->_f_cnt          = tm._f_cnt;

                this->_v_pos          = tm._v_pos;
                this->_v_n            = tm._v_n;
                this->_e_len          = tm._e_len;
                this->_e_len_average  = tm._e_len_average;
                this->_f_n            = tm._f_n;
                this->_f_center       = tm._f_center;
                this->_f_size         = tm._f_size;

                this->_center                   = tm._center;
                this->_bounding_box[0]          = tm._bounding_box[0];
                this->_bounding_box[1]          = tm._bounding_box[1];
                this->_bounding_box_diagonal    = tm._bounding_box_diagonal;
                this->_bounding_sphere_radius   = tm._bounding_sphere_radius;

                return *this;
            }

            void Init(void)
            {
                _v_2_e.clear();
                _v_2_f.clear();
                _e_2_v.clear();
                _e_2_f.clear();
                _f_2_v.clear();
                _f_2_e.clear();

                _v_cnt = 0;
                _e_cnt = 0;
                _f_cnt = 0;

                _v_pos.clear();
                _v_n.clear();
                _e_len.clear();
                _e_len_average = 0.0f;
                _f_n.clear();
                _f_center.clear();
                _f_size.clear();

                _center = math::vec3<T>::Zero();
                _bounding_box[0] = math::vec3<T>::Zero();
                _bounding_box[1] = math::vec3<T>::Zero();
                _bounding_box_diagonal = 0.0f;
                _bounding_sphere_radius = 0.0f;
            }

            //------------------------------------------------------------------------------------------------------------
            // Load/Save 
            //------------------------------------------------------------------------------------------------------------
            bool IsLoaded(void)
            {
                return _v_cnt;
            }

            bool LoadFromFile_OBJ(const std::string& filename)
            {
                this->Init();

                FILE *fp = fopen(filename.c_str(),"r");
                if (fp == NULL) {
                    return false;
                }

                char line[1024];
                char buf[1024];
                T d[3];
                std::vector<int> i(3);

                while (fgets(line, 1024, fp)) {
                    buf[0] = 0;
                    sscanf(line, "%s", buf);

                    if (std::string(buf) == std::string("v")) {
                        char fmt[1024];
                        const char* p = math::get_type_indicator<T>();
                        sprintf(fmt, "%%s %%%s %%%s %%%s", p, p, p);
                        sscanf(line, fmt, buf, &d[0], &d[1], &d[2]);
                        this->_v_pos.push_back(math::vec3<T>(d));
                    } else if (std::string(buf) == std::string("f")) {
                        std::vector<int> t(3);
                        std::vector<int> tt(3);
                        char *p = strchr(line, '/');

                        if (p == NULL) {
                            sscanf(line, "%s %d %d %d", buf, &i[0], &i[1], &i[2]);
                        } else if (*(++p) == '/') {
                            sscanf(line, "%s %d//%d %d//%d %d//%d", buf, &i[0], &t[0], &i[1], &t[1], &i[2], &t[2]);
                        } else {
                            sscanf(line, "%s %d/%d/%d %d/%d/%d %d/%d/%d", buf, &i[0], &t[0], &tt[0], &i[1], &t[1], &tt[1], &i[2], &t[2], &tt[2]);
                        }
                        for (int j = 0; j < 3; j++) {
                            if (i[j] < 0) {
                                i[j] = -i[j] - 1;
                            } else {
                                i[j]--;
                            }
                        }
                        this->_f_2_v.push_back(i);
                    }
                }
                fclose(fp);

                _v_cnt = _v_pos.size();
                _f_cnt = _f_2_v.size();

                this->CalcuTopo();
                this->CalcuGeometry();

                return true;
            }

            bool SaveToFile_OBJ(const std::string& filename)
            {
                FILE *fp = fopen(filename.c_str(),"w");
                if (fp == NULL) {
                    return false;
                }

                // save vertex position
                char fmt[1024];
                const char* p = get_type_indicator<T>();
                sprintf(fmt, "v %%.8%s %%.8%s %%.8%s\n", p, p, p);
                for (int i = 0; i < _v_cnt; i++) {
                    fprintf(fp, fmt, _v_pos[i][0], _v_pos[i][1], _v_pos[i][2]);
                }

                // save face info
                for (int i = 0; i < _f_cnt; i++) {
                    fprintf(fp, "f");
                    for (int j = 0; j < _f_2_v[i].size(); j++) {
                        fprintf(fp, " %d", _f_2_v[i][j] + 1);
                    }
                    fprintf(fp, "\n");
                }
                fclose(fp);

                return true;
            }

            bool LoadFromFile_TM(const std::string& filename)
            {
                this->Init();

                FILE *fp = fopen(filename.c_str(),"r");
                if (fp == NULL) {
                    return false;
                }

                int zz, i;
                fscanf(fp,"%d %d %d",&_v_cnt,&_f_cnt,&zz);
                _v_pos.resize(_v_cnt);
                _f_2_v.resize(_f_cnt);
                // load vertex position
                char fmt[1024];
                const char* p = get_type_indicator<T>();
                sprintf(fmt, "%%%s %%%s %%%s", p, p, p);
                for (i=0;i<_v_cnt;++i)
                    fscanf(fp, fmt, &_v_pos[i].x, &_v_pos[i].y, &_v_pos[i].z);
                // load face info
                for (i=0;i<_f_cnt;++i) {
                    _f_2_v[i].resize(3);
                    fscanf(fp,"%d %d %d",&_f_2_v[i][0],&_f_2_v[i][1],&_f_2_v[i][2]);
                }
                fclose(fp);

                this->CalcuTopo();
                this->CalcuGeometry();

                return true;
            }

            bool SaveToFile_TM(const std::string& filename)
            {
                FILE *fp = fopen(filename.c_str(), "w");
                if (fp == NULL) {
                    return false;
                }

                fprintf(fp, "%d %d %d\n", _v_cnt, _f_cnt, _f_cnt);
                // save vertex position
                char fmt[1024];
                const char* p = get_type_indicator<T>();
                sprintf(fmt, "%%.8%s %%.8%s %%.8%s\n", p, p, p);
                for (int i = 0; i < _v_cnt; ++i)
                    fprintf(fp, fmt, _v_pos[i].x, _v_pos[i].y, _v_pos[i].z);
                // save face info
                for (int i = 0; i < _f_cnt; ++i) {
                    fprintf(fp, "%d %d %d\n", _f_2_v[i][0], _f_2_v[i][1], _f_2_v[i][2]);
                }
                fclose(fp);

                return true;
            }

            //------------------------------------------------------------------------------------------------------------
            // Topology Informations
            //------------------------------------------------------------------------------------------------------------
            std::vector<std::vector<int> > _v_2_e;
            std::vector<std::vector<int> > _v_2_f;
            std::vector<std::vector<int> > _e_2_v;
            std::vector<std::vector<int> > _e_2_f;
            std::vector<std::vector<int> > _f_2_v;
            std::vector<std::vector<int> > _f_2_e;

            int _v_cnt;
            int _e_cnt;
            int _f_cnt;

            void CalcuTopo()
            {
                _v_cnt = (int) _v_pos.size();
                _f_cnt = (int) _f_2_v.size();
                _e_cnt = 0;

                int i, j, k;
                int l1, l2, l3, p1, p2, p3;

                _rprt_v_non_edge.clear();
                _rprt_v_non_face.clear();
                _rprt_e_non_face.clear();
                _rprt_e_multi_face.clear();

                // <--
                // _v_2_f
                _v_2_f = std::vector<std::vector<int> >(_v_cnt);

                for (i = 0; i < _f_cnt; ++i) {
                    for (j = 0; j < 3; ++j)
                        _v_2_f[_f_2_v[i][j]].push_back(i);
                }

                // <--
                // v_2_v
                std::vector<std::vector<int> > v_2_v = std::vector<std::vector<int> >(_v_cnt);

                for (i = 0; i < _f_cnt; ++i) {
                    p1 = _f_2_v[i][0];
                    p2 = _f_2_v[i][1];
                    p3 = _f_2_v[i][2];
                    v_2_v[p1].push_back(p2);
                    v_2_v[p1].push_back(p3);
                    v_2_v[p2].push_back(p3);
                    v_2_v[p2].push_back(p1);
                    v_2_v[p3].push_back(p1);
                    v_2_v[p3].push_back(p2);
                }
                for (i = 0; i < _v_cnt; ++i) {
                    // validation
                    assert(v_2_v[i].size() % 2 != 1);
                    if (v_2_v[i].size() == 1) {
                        _rprt_v_non_face.push_back(i);
                        continue;
                    } else if (v_2_v[i].size() == 0) {
                        _rprt_v_non_edge.push_back(i);
                        continue;
                    }

                    const std::vector<int>::iterator begin = v_2_v[i].begin(), end = v_2_v[i].end();
                    std::sort(begin, end);
                    v_2_v[i].erase(std::unique(begin, end), end);
                }
                // validation
                for (i = 0; i < _v_cnt; ++i) {
                    for (j = 0; j < (int) v_2_v[i].size(); ++j) {
                        assert(v_2_v[i][j] != i); // a vertex should never connect to itself
                    }
                }

                // <--
                // _v_2_e
                // _e_2_v
                _v_2_e = std::vector<std::vector<int> >(_v_cnt);
                for (i = 0; i < _v_cnt; ++i) {
                    _v_2_e[i] = std::vector<int>(v_2_v[i].size(), -1);
                }
                _e_2_v.clear();

                {
                    std::vector<int> e(2);
                    for (i = 0; i < _v_cnt; ++i) {
                        e[0] = p1 = i;
                        for (j = 0; j < (int) v_2_v[i].size(); ++j) {
                            e[1] = p2 = v_2_v[p1][j];
                            if (p1 > p2) {
                                continue;
                            }
                            // add _e_2_v
                            _e_2_v.push_back(e);
                            // add _v_2_e
                            _v_2_e[p1][j] = _e_cnt;
                            for (k = 0; k < (int) v_2_v[p2].size(); ++k) {
                                if (v_2_v[p2][k] == p1) {
                                    _v_2_e[p2][k] = _e_cnt;
                                    break;
                                }
                            }
                            assert(k != v_2_v[p2].size()); // p2 should have p1 in its v_2_v
                            _e_cnt++;
                        }
                    }
                }
                // validation
                for (i = 0; i < _e_cnt; i++) {
                    assert(_e_2_v[i].size() == 2); // each edge should contain two vertices
                }
                for (i = 0; i < _v_cnt; i++) {
                    for (j = 0; j < (int) _v_2_e[i].size(); ++j) {
                        assert(_v_2_e[i][j] != -1); // _v_2_e should be filled
                    }
                }

                // <--
                // _f_2_e
                // _e_2_f
                _f_2_e = std::vector<std::vector<int> >(_f_cnt);
                _e_2_f = std::vector<std::vector<int> >(_e_cnt);

                for (i = 0; i < _f_cnt; ++i) {
                    for (l1 = 0; l1 < 3; ++l1) {
                        l2 = (l1 + 1) % 3;
                        l3 = (l1 + 2) % 3;
                        p1 = _f_2_v[i][l1];
                        p2 = _f_2_v[i][l2];
                        p3 = _f_2_v[i][l3];

                        for (j = 0; j < (int) v_2_v[p1].size(); ++j) {
                            if (v_2_v[p1][j] == p2) {
                                break;
                            }
                        }
                        assert(j != v_2_v[p1].size()); // vertices in a face should be connected

                        int eindex = _v_2_e[p1][j];
                        _f_2_e[i].push_back(eindex);
                        _e_2_f[eindex].push_back(i);
                    }
                }

            }

            //------------------------------------------------------------------------------------------------------------
            // Geometry Informations
            //------------------------------------------------------------------------------------------------------------
            std::vector<math::vec3<T> >  _v_pos;
            std::vector<math::vec3<T> >  _v_n;
            std::vector<T>              _e_len;
            T                           _e_len_average;
            std::vector<math::vec3<T> >  _f_n;
            std::vector<math::vec3<T> >  _f_center;
            std::vector<T>              _f_size;

            math::vec3<T>                _center;
            math::vec3<T>                _bounding_box[2];
            T                           _bounding_box_diagonal;
            T                           _bounding_sphere_radius;

            void CalcuGeometry()
            {
                int i, j;

                // _f_n[], _f_size[], f_c[]
                _f_n.resize(_f_cnt);
                _f_size.resize(_f_cnt);
                _f_center.resize(_f_cnt);
                for (i=0;i<_f_cnt;++i) {
                    math::vec3<T> v0 = _v_pos[_f_2_v[i][0]];
                    math::vec3<T> v1 = _v_pos[_f_2_v[i][1]];
                    math::vec3<T> v2 = _v_pos[_f_2_v[i][2]];
                    math::vec3<T> zz = (v1 - v0) ^ (v2 - v0);
                    if (zz == math::vec3<T>::Zero()) {
                        _f_n[i] = math::vec3<T>::Zero();
                        _f_size[i] = 0.0f;
                    } else {
                        _f_size[i] = zz._norm_ex(_f_n[i]) / 2;
                    }
                    _f_center[i] = (v0 + v1 + v2) * (one<T>() / 3);
                }

                // _v_n[]
                _v_n.resize(_v_cnt);
                for (i = 0; i < _v_cnt; ++i) {
                    _v_n[i] = math::vec3<T>::Zero();
                    for (j=0;j<(int) _v_2_f[i].size();++j) {
                        int findex = _v_2_f[i][j];
                        _v_n[i] += _f_n[findex] * _f_size[findex];
                    }
                    if (_v_n[i] != math::vec3<T>::Zero()) {
                        _v_n[i]._norm();
                    }
                }

                // _bounding_box[], _bounding_box_diagonal
                _bounding_box[0] = math::vec3<T>(math::inf<T>(), math::inf<T>(), math::inf<T>());
                _bounding_box[1] = math::vec3<T>(-math::inf<T>(), -math::inf<T>(), -math::inf<T>());
                for (i = 0; i < _v_cnt; ++i) {
                    for (j = 0; j < 3; j++) {
                        _bounding_box[0][j] = math::min(_bounding_box[0][j], _v_pos[i][j]);
                    }
                    for (j = 0; j < 3; j++) {
                        _bounding_box[1][j] = math::max(_bounding_box[1][j], _v_pos[i][j]);
                    }
                }
                _bounding_box_diagonal = (_bounding_box[0] - _bounding_box[1]);

                // _center
                _center = (_bounding_box[0] + _bounding_box[1]) * (one<T>() / 2);

                // _bounding_sphere_radius
                _bounding_sphere_radius = 0.0f;
                for (i = 0; i < _v_cnt; ++i) {
                    T d = _center - _v_pos[i];
                    _bounding_sphere_radius = math::max(_bounding_sphere_radius, d);
                }

                // _e_len[], _e_len_average
                _e_len.resize(_e_cnt);
                _e_len_average = 0;
                for (i = 0;i<_e_cnt;++i) {
                    _e_len_average += (_e_len[i] = (T) (_v_pos[_e_2_v[i][0]] - _v_pos[_e_2_v[i][1]]));
                }
                _e_len_average /= _e_cnt;

            }

            T Error(const tmesh& another) const
            {
                assert(_v_cnt == another._v_cnt);

                T err = 0.0f;

                for (int i = 0; i < _v_cnt; i++) {
                    err += (T) (_v_pos[i] - another._v_pos[i]);
                }

                return err / _v_cnt;
            }

            T RelativeError(const tmesh& another) const
            {
                return Error(another) / _bounding_sphere_radius;
            }


            //------------------------------------------------------------------------------------------------------------
            // Debug 
            //------------------------------------------------------------------------------------------------------------
            friend std::ostream& operator << (std::ostream& os, tmesh& tm)
            {
                os << "[ tmesh" << std::endl;
                os << "    " << (int) tm._v_cnt << " vertices" << std::endl;
                os << "    " << (int) tm._f_cnt << " faces" << std::endl;
                os << "    " << (int) tm._e_cnt << " edges" << std::endl;
                os << "]" << std::endl;

                return os;
            }

        public:
            //------------------------------------------------------------------------------------------------------------
            // Report Informations
            //------------------------------------------------------------------------------------------------------------
            std::vector<int> _rprt_v_non_edge;        // an isolated vertex
            std::vector<int> _rprt_v_non_face;        // with only one edge
            std::vector<int> _rprt_e_non_face;        // with no face
            std::vector<int> _rprt_e_multi_face;      // with multiple faces
        };

    }

}
